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Abstract 
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Due to their structural flexibility, spacecraft and space manipulators are multibody 
systems with complex dynamics and possess a large number of degrees of freedom. This 
publication uses the spatial operator algebra methodology to develop a new dynamics for- 
mulation and spatially recursive algorithms for such flexible multibody systems. A key 
feature of the formulation is that the operator description of the flexible system dynamics is 
identical in form to the corresponding operator description of the dynamics of rigid multi- 
body systems. A significant advantage of this unifying approach is that it allows ideas and 
techniques for rigid multibody systems to be easily applied to flexible multibody systems. 
The algorithms use standard finite-element and assumed modes models for the individual 
body deformation. 

A Newton-Euler Operator Factorization of the mass matrix of the multibody system 
is first developed. It forms the basis for recursive algorithms such as for the inverse dynam- 
ics, the computation of the mass matrix, and the composite body forward dynamics for the 
system. Subsequently, an alternative Innovations Operator Factorization of the mass matrix, 
each of whose factors is invertible, is developed. It leads to an operator expression for the 
inverse of the mass matrix, and forms the basis for the recursive articulated body forward 
dynamics algorithm for the flexible multibody system. For simplicity, most of the develop- 
ment in this publication focuses on serial chain multibody systems. However, extensions of 
the algorithms to general topology flexible multibody systems are described. 

While the computational cost of the algorithms depends on factors such as the topol- 
ogy and the amount of flexibility in the multibody system, in general, it appears that in 
contrast with the rigid multibody case, the articulated body forward dynamics algorithm is 
the more efficient algorithm for flexible multibody systems containing even a small number 
of flexible bodies. The variety of algorithms described here permits a user to choose the al- 
gorithm which is optimal for the multibody system at hand. The availability of a number 
of algorithms is even more important for real-time applications, where implementation on 
parallel processors or custom computing hardware is often necessary to maximize speed. 
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1. Introduction 


Due to their structural flexibility, spacecraft and space manipulators are multibody systems 
with complex dynamics and possessing a large number of degrees of freedom (dofs). This 
publication describes a spatial operator formulation for the analysis and development of 
efficient dynamics algorithms for such flexible multibody systems. This approach represents 
a reformulation and extension of the spatial operator algebra methodology [1] for the dynamics 
of flexible multibody systems [2]. 

The formulation developed in this publication uses spatial operators and closely par- 
allels the corresponding formulation for rigid multibody systems. In effect it represents a 
unified methodology for the analysis and development of efficient dynamics algorithms for 
general topology, rigid/flexible multibody systems. Indeed, the analysis and recursive algo- 
rithms developed here for the inverse dynamics, the computation of the mass matrix, and 
the forward dynamics of the flexible multibody system closely resemble the corresponding 
analysis and algorithms for the rigid multibody system. 

It is assumed that all deformations of the bodies are small and a linear model of 
elasticity is used. However, large articulation at the hinges is allowed. No special assumptions 
are made regarding the nature of the component bodies, and the algorithms use standard 
finite-element and assumed modes models for the body flexibility. For notational simplicity, 
and without any loss in generality, the main focus of this publication is on flexible multibody 
systems with serial chain topology. The extensions required for systems with tree and closed- 
chain topology are discussed at the end of the publication. 

Section 2 summarizes the symbols and notation used in this publication. Section 
3 contains the development of the equations of motion for the multibody system. The 
recursive relationships between the modal spatial velocities, modal spatial accelerations, and 
modal spatial forces are described. Using these, spatial operators are used to develop the 
Newton-Euler Operator Factorization of the system mass matrix. 

Section 4 describes a recursive Newton-Euler inverse dynamics algorithm for the 
flexible multibody system. This algorithm computes the vector of generalized forces cor- 
responding to a given state and the vector of generalized accelerations for the multibody 
system. In Section 5, the Newton-Euler Operator Factorization of the mass matrix is used 
to develop the recursive composite body forward dynamics algorithm for the system. A part 
of this algorithm consists of an algorithm for the recursive computation of the multibody 
system mass matrix. 

Section 6 describes operator factorization and inversion results that form the basis for 
the recursive articulated body forward dynamics algorithm. First, a recursive algorithm for 
the computation of certain articulated body quantities is defined. These quantities are used 
to develop a new operator factorization, denoted the Innovations Operator Factorization of 
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the system mass matrix. In contrast with the Newton-Euler Operator Factorization, the 
factors in the Innovations Operator Factorization are square and invertible operators. This 
fact is used to develop an operator expression for the inverse of the mass matrix. A more 
efficient algorithm for computing the articulated body quantities is described at the end of 
the section. Using the operator expression for the mass matrix inverse, Section 7 describes 
the recursive articulated body forward dynamics algorithm for the multibody system. This 
algorithm requires neither the explicit formation of the system mass matrix nor its inversion. 

Section 8 discusses the computational complexity of the various algorithms described 
in this publication. It is shown that the articulated body forward dynamics algorithm is 
more efficient than the composite body forward dynamics algorithm, even for systems with 
a modest number of bodies and flexible modes. Section 9 describes the extensions of the 
dynamics formulation and algorithms developed here to flexible multibody systems with tree 
and closed-chain topologies. 


2. Notation 


Coordinate free spatial notation [1, 3] is used throughout this publication. Briefly a spatial 
velocity of a frame is a 6-dirnensional quantity, whose upper 3 elements are the angular 
velocity while the lower 3 elements are the linear velocity of the frame. A spatial force is a 
6-dimensional quantity, whose upper 3 elements correspond to a moment vector while the 
lower 3, to a force vector. 

In the stacked notation used in this publication, a variety of indices are used to 
identify different spatial quantities. Thus for instance, V 3 (jk) denotes the spatial velocity of 
the j th node on the k th body, V s (k ) = col|V*(jfc)j denotes the vector of the spatial velocities 

of all the nodes on the k th body, while V s = coljl^(A;)j denotes the vector of the spatial 
velocities of all the nodes for all the bodies in the serial chain. The index k will be used to 
refer to both the k th body as well as the k th body reference frame with the usage being 
apparent from the context. Some key quantities used in this publication are defined below 
(also see Figure 1). 
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Figure 1: Illustration of links and hinges in a flexible serial multibody system 


General Quantities: 

z — [z] x 6 7£ 3x3 and denotes the skew-symmetric cross-product matrix associ- 
ated with the 3-dimensional vector z 

dx 

x = — and denotes the time derivative of x with respect to an inertial frame 
d t 

x denotes the time derivative x with respect to the body-fixed (rotating) frame 


diag jx(Ar)J 
col{x(fc)j 

K x ,y) 


denotes a block diagonal matrix whose k th diagonal element is x(k) 
denotes a column vector whose k th element is x(k) 

£ H 3 , the vector from point x to point y 




I l(x,y) 


£ 7Z 6 * 6 , the spatial transformation operator which 


V 0 / I . . 

forms spatial velocities and forces between points x and y 

Individual Body Nodal Data: 

n 3 (k ) number of nodes on the k th body 

Tk the body reference frame with respect to which the deformation field for the 
k th body is measured. The motion of this frame characterizes the motion of 
the k th body as a rigid body. 


jk denotes the j th node on the k th body 

^o(^-ijk) G 7£ 3 , the vector from Tk to the location (before deformation) of the j th node 
reference frame on the k th body 

Si(j k ) € H 3 , the translational deformation of the j th node on the k th body 

— lo(kijk) + $l(jk) G the vector from Tk to the location (after deformation) 
of the j th node reference frame on the k th body 


^w(ifc) G ^ 3 7 the deformation angular velocity of the j th node on the k th body with 
respect to the body frame Tk 

S v (jk) G 7^ 3 , the deformation linear velocity of the j th node on the k th body with 
respect to the body frame Tk 


u(jk) 


£ 7Z e , the spatial displacement of node jk . The translational component of 
u(jk) is Si(jk ), while its time derivative with respect to the body frame JF*. is 


u(jk) 


1 Uh ) ' 

^ f>v{jk) j 
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JUk) 
p{. h) 
™(jk) 
M.Vk) 

M,(k) 

K,(k) 


G 72. 3x3 , the inertia tensor about the nodal reference frame for the j th node 
on the k th body 

(E Tl 3 , the vector from the nodal reference frame to the node center of mass 
for the j th node on the k th body 

the mass of the j th node on the k th body 


J(jk) m(j k )p(j k ) ' 


V 


€ 7?. 6x6 , the spatial inertia about the nodal 

~™(jk)p(jk) rn(j k )I ) 
reference frame for the j th node on the k th body 

diag{M 4 (ifc)} G ^*n»(*)x6n,(fc)^ structural mass matrix for the k th body 
e 'J^6n s (k)x6n,(k)^ s t ruc tural stiffness matrix for the k th body 


Individual Body Modal Data: 


n m {k ) 
~Af(k) 

ri(k) 

a m 


the number of assumed modes for the k th body 

= n m (k) + 6, the number of deformation plus rigid-body dofs for the k th body 

£ %^m(k)^ ^g vec tor of modal deformation variables for the k th body 

G Tl 3 , the modal slope (or differential change in orientation) displacement 
vector for the r th mode at the j l k h nodal reference frame. 


\i(k) — G 7£ 3xn "d fc ), the modal slope displacement influence 

vector for the j th node due to all the modes for the k ih body. Note that 

Mi*) = V(k)f){k). 

7 J(j k) e 'll 3 , the modal translational displacement vector for the r ih mode at the j[ h 

nodal reference frame 


l 3 (k) 


n m 


ip(jfc) 


= [7i(^)> " ' >7n (*)(^)1 € 7£ 3xnm M the modal translational displacement in- 
fluence vector for the j ttl node due to all the modes for the k tfl body. Note 
that 6,{j k ) = l j (k)ri{k) and 5„(i*) = ~f 3 (k)fi(k). 


/ H(k) X 


V 


G 'll 6 , the modal spatial displacement vector for the r th mode 

l 3 r{k) r 

at the j[ h nodal reference frame 


= [ni«, •••,<.,*,(*) 

the j[ h node. The spatia 


G 7£ 6xn ">W 

deformation 


, the modal spatial influence vector for 
of node jk is given by u(jk) = bP (k)r](k). 
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IT(fc) = coljrP(fc)j G '!Z 6n, W* nm ( k \ the modal matrix for the k ih body. The r th 
column of II(fc) is denoted II r (k) G 7Z 6n ’^ and is the mode shape function for 
the r th assumed mode for the k th body. The deformation field for the k th body 
is given by u(k ) = n(fc)7/(fc), while u(k ) = I1(A:)^(A:). 

M m {k) G 7 ^ mo ^ a i mass matrix for the k th body. 

Km(k) G the modal stiffness matrix for the k th body. 

Multibody Data: 


N 

jr 


n r {k) 

AT(k) 

X 

d k 

t k 

O k 


Ot 

m 

m 

A v{k) 


H'(k) 

m 


number of bodies in the serial flexible multibody system 

= XjtLi A (A:), the overall dofs in the serial chain obtained by disregarding the 
hinge constraints 

the number of dofs for the k th hinge 

= n m (k ) + n T (k), the number of deformation plus hinge dofs for the k th body 

= HfcLi ■N’(k), the overall deformation plus hinge dofs for the serial chain 

denotes the node on the k th body to which the k th hinge is attached 

denotes the node on the k th body to which the (k — l)*' 1 hinge is attached 

the reference frame for the k th hinge on the k th body. This frame is fixed to 
node d k . 

the reference frame for the k th hinge on the (k - f 1)"* body. This frame is fixed 
to node fjt+i- 

G 7 Z nr ( k \ the vector of configuration variables for the k th hinge 
G TZ n ^ k \ the vector of generalized speeds for the k th hinge 


1 A„(l') ) 


G 7l e , the relative spatial velocity for the k th hinge defined as 


\ A .(A) j 

the spatial velocity of frame O k with respect to frame 0~l 

G 7£ 6xn d*) 5 the joint map matrix for the k th hinge such that Av(&) = H*(k)fl(k). 


( rj(k) 

G 'Rr' k \ the vector of (deformation plus hinge) generalized con- 

*(*) / 

figuration variables for the k th body 
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m 


£ the vector of (deformation plus hinge) generalized 


x(k) 


V(k) 


V(O t ) 

V(Ot) 

Vs(jk) 

V m (k) 


^m(^) 

a m (k) 

b m (k) 

f m (k) 

fs(jk) 

m 

T( k) 
H?(k) 

H{k) 


( 


\ 


m ) 

motion variables (or generalized speeds) for the k th body 


= v{r„) = 


u }(k) 


6 R 6 , the spatial velocity of the k th body reference 


\ v ( k ) j 

frame Tk, with w(Jfc) and v(k ) denoting the angular and linear velocities re- 
spectively of frame Tk 

G R 6 , the spatial velocity of frame Ok 

G R 6 , the spatial velocity of frame Of 

G It 6 , the spatial velocity of the j th node on the k th body. 

G R 6 , the spatial acceleration of the j th node on the k th body. 


/ m x 


V(k) 


) 


G RF( k \ the modal spatial velocity of the k th body 


= K„(Jt) G RF^ k \ the modal spatial acceleration of the k th body 

G RF^ k \ the modal Coriolis and centrifugal accelerations for the k th body 

G RF( k \ the modal gyroscopic forces for the k th body 

G RF^ k \ the modal spatial force of interaction between the k th and 
( k + l) th bodies 

G R 6 , the spatial force at node jk 
G R 6 , the effective spatial force at frame Tk 
G R^ k \ the generalized force for the k th body 

= H(k)(f>(Ok,k) G 72” r (* : ) x6 , the joint map matrix referred to frame Tk for the 
k th hinge 


' i -tn'(t)]- x 


'j^N’(k)y.U(k) ^ (Reformation plus hinge) modal joint 


\,0 Hr(k) J 

map matrix for the k th body 
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A(k) 


G 7^(*)x« j relates spatial forces and velocities between node 


B(k + l,k) 


*(k+l,k) 


= ' tn’(*)J* ' 

v <KM*) j 

tk and frame Tk 


= [0, <f>(tk+i,k)] G 72. 6x ^( fc ), relates spatial forces and velocities between node 
tjt+i <md frame Tk 


= ,4(fc + !)*(*+ 1,*) = 


( o [n‘(* + i)]V(4+i,fc) \ 


y 0 <£(& + 1, A;) J 

interbody transformation operator which relates modal 
locities between the k th and ( k + l) th bodies 


6 the 

spatial forces and ve- 


( 


0 


\ 


C(k,k~ 1 ) 


4>{tk, k — 1) 


g 7?, 6n *W x6 


\ 


0 


/ 


5(ib) = It), <j>(k,2k),- ■ ■ , 4>(k, n,(fc))] G 72. 6x6 "*(*), relates the spatial velocity 

of frame Tk to the spatial velocities of all the nodes on the k th body when the 
body is regarded as being rigid 

M € ft^ xAr , the multibody system mass matrix 

C G , the vector of Coriolis, centrifugal and elastic forces for the multibody 
system 


3. Equations of Motion for Flexible Serial Chains 


In this section we develop the equations of motion for a serial flexible multibody system 
consisting of N flexible bodies. Recursive relationships between the modal spatial velocities, 
accelerations and forces are developed. Spatial operators are introduced to express these 
relationships in a compact form, and obtain the Newton-Euler Operator factorization of the 
mass matrix for the multibody system. 

It is assumed that each flexible body has a lumped mass model with a a rigid body 
being located at each node. The number of nodes on the k th body is denoted n a (k). The 
j th node on the k th body is referred to as the jj. h node. Each body has associated with it 
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a body reference frame, denoted Tk for the k tfl body. The deformation of the nodes on the 
body is described with respect to this body reference frame, while the rigid body motion of 
the k th body is characterized by the motion of frame Tk- 

The 6-dimensional spatial deformation (slope plus translational) of node jk (with 
respect to frame f k ) is denoted u(j k ) E ft 6 . The overall deformation field for the k th body is 
defined as the vector u[k) = col{u(jjt)| E 7^ 6n d fc ). The vector from frame iF k to the reference 
frame on node jk is denoted l(k,jk ) E ft 3 - 

The spatial inertia of the j th node is defined as 




J(h) 

K - m{j k )P{jk ) 


™{jk)p{jk) ^ 
™{jk)I j 


e ft 6 * 6 


(3.1) 


where ffi^Ji t), PO'fci) <md m(jk') are the inertia tensor about the node reference frame, the 
vector from the node reference frame to its center of mass, and the mass, respectively, for 
thej tft node on the k th body. The structural mass matrix for the k th body M„{k ) is the block 
diagonal matrix 

M,(k) = diag{M a (i,)} € 7l 6 ”‘ (fc)x6Mfc) (3-2) 

The structural stiffness matrix is denoted K a (k ) E 7£ 6n 4 fc ) x6n *( <: ). 

As shown in Figure 1, the bodies in the serial chain are numbered in increasing 
order from tip to base. We use the terminology inboard (outboard) to denote the direction 
along the serial chain towards (away from) the base body. The k th body is attached on 
the inboard side to the (k + l) th body via the k th hinge, and on the outboard side to the 
( k - l) th body via the ( k - l) t/l hinge. On the k th body, the node to which the outboard 
hinge (the ( k - l) th hinge) is attached is referred to as node t k , while the node to which the 
inboard hinge (the k th hinge) is attached is denoted node d k . Thus the k ih hinge couples 
together nodes d k and t k+ Attached to each of these nodes are the k th hinge reference frames 
denoted O k and O k , respectively. The number of dofs for the k th hinge is denoted n r {k). The 
vector of configuration variables for the k th hinge is denoted &(k) E 'R- nr ^\ while its vector 
of generalized speeds is denoted /3(k) E In general, when there are nonholonomic 

hinge constraints, the dimensionality of f3(k) may be less than that of 6(k). For notational 
convenience, and without any loss in generality, we assume here that the dimensions of the 
vectors 8(k) and (3(k) are equal. In most situations /3(k) is simply 0. However there are 
many cases where the use of quasicoordinates simplifies the dynamical equations of motion 
and an alternative choice for (3(k ) may be preferable. The relative spatial velocity A v{k) 
across the hinge is given by H*(k){3(k), where H*(k) denotes the joint map matrix for the 
k th hinge. 
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We assume that there exists a set of n m (k ) assumed modes for the k th body, and let 
n^(Jfc) £ 7Z 6 denote the modal spatial displacement vector at the j k h node for the r th mode. 
We also define the modal spatial displacement influence vector IP (k) £ 7 ^exn m (fc) f or jth 
node and the modal matrix II(£) £ 7 ^6n»(t)xn m (*:) f or j^th jj 0( jy as follows: 

nJ(i)=[ni(t),.-,ni. w (t)] and n(t) = a>i{ip(t)} 

The r th column of Ft (k) is denoted U r (k) and defines the mode shape for the r th assumed 
mode for the k th body. With 7}(k ) € 7 Z nrn ^ denoting the vector of modal deformation 
variables for the k th body, the spatial deformation of node j \ t and the spatial deformation 
field u(k) for the k th body are given by 

u(j k ) = U j (k)7](k) and u{k) = U(k)rj(k) (3.3) 


In the multibody context, it is often convenient to choose the k th body reference frame 
T k as fixed to node d k at the inboard hinge. For this choice, the assumed modes are referred 
to as cantilever modes, and for which 

(A;) = 0 and r = 1 - - • n m (k) (3.4) 


As a consequence, node d k exhibits zero deformation (u(d k ) = 0). Alternative choices of 
modes, which are often preferred for control analysis and design, are the free-free modes . For 
this case, the reference frame T k is not fixed to any node, but is rather assumed to be fixed 
to the undeformed body, so that all nodes exhibit nonzero deformation. We assume here 
that when free-free assumed modes are used, only the deformation modes, and none of the 
rigid body modes, have been included in the mode set. The dynamics model and algorithms 
developed here handle both types of modes’ cases, with some additional computational 
simplifications arising from Eq. (3.4) when cantilever modes are used. For a related discussion 
regarding the choice of reference frame and modal representations for a flexible body see 
reference [4]. 


The vector of generalized configuration variables d(k) and generalized speeds x(^) f° r 
the k th body are defined as 


i 9(k) = 


m j 


£ 7l^ k) and x(jfc) = 


' m ^ 
\ m j 


£ n" ik) 


(3.5) 


where Af(k) = n m (k ) + n r (k). The overall vectors of generalized configuration variables d 
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and generalized speeds \ for the serial multibody system are defined as 

i? = col{tf(Jk)} € and X = col{x(*)} 6 K* (3.6) 

where Af = X2JfcLi Af{k). The number of overall dofs for the multibody system is Af . The 
state of the multibody system is defined by the pair of vectors {tf, x}- For a given system 
state {$, x}, the equations of motion define the relationship between the vector of generalized 
accelerations x and the vector of generalized forces T G for the system. The inverse 
dynamics problem consists of computing the vector of generalized forces T for a prescribed 
set of generalized accelerations The forward dynamics problem is the converse one and 
consists of computing the set of generalized accelerations x resulting from a set of generalized 
forces T. The equations of motion for the system are developed in the remainder of this 
section. In Section 4, a recursive Newton-Euler inverse dynamics algorithm is developed, 
while Section 5 and Section 7 describe two alternative forward dynamics algorithms. 


3.1 Recursive Propagation of Velocities 


In this section, recursive relationships as well as operator expressions for the modal spatial 
velocities of the bodies in the serial chain are developed. 

Let V(k) denote the spatial velocity of the k th body reference frame i.e., 


V(k) = 


u>(k) 

K v(k) j 


<e n 6 


where u>(k) and v(k) denote the angular and linear velocities respectively of jF/t- The spatial 
velocity G 77 6 of node tk+ 1 (on the inboard of the k th hinge) is related to the spatial 

velocity V(k + 1) of the ( k + l) i/l body reference frame .Fjt+i, and the modal deformation 
variable rates i]{k + 1) as follows: 


T»(Lt+i) — + lj^fc+i)l r (^ + 1) + w(<fc + i) 

= f*(k + Mfc+i )V{k + 1) + n l (k + 1 )y(k + 1) 


The spatial transformation operator 4>(x,y) G 77. 6x6 above is defined to be 


4>( x ,y) 


f I l{x,y) ^ 
K° 1 ) 


(3.7) 


(3.8) 
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where l(x,y) € 7Z 3 denotes the vector between the points x and y. Note that the following 
important (group) property holds: 


<t>(x,y)<f>(y,z) = <j>(x,z) 

for arbitrary points x, y and As in Eq. (3.7), and all through this publication, the index 
k will be used to refer to both the k th body as well as to the k th body reference frame T k 
with the specific usage being evident from the context. Thus V(k ) may be read as V{T\ t), 
and 4>{k,t k ) as <f>(Tk,t k ). 

The spatial velocity V(O k ) of frame O k on the inboard side of the k th hinge is related 
to V a (t k+ 1 ) via 


V(Ot) = f(t k+1 ,O k )V s (t k+1 ) (3.9) 

The relative spatial velocity A v(k) across the k th hinge is given by H*(k)(3(k), and so the 
spatial velocity V(O k ) of frame O k on the outboard side of the k th hinge is 

V(O k ) = V(Ot) + H*(k)/3(k ) (3.10) 

The spatial velocity V(k) of the k th body reference frame is given by 

V(k) = <f>*(O k ,k)V(O k ) - u(d k ) = <f>*{O k ,k)V{O k ) - n d (k)r,(k) (3.11) 

Putting together Eq. (3.7), Eq. (3.9), Eq. (3.10) and Eq. (3.11), it follows that 

V(k) = <j>*(k + l,k)V(k+l) + <j>*(t k+1 ,k)W{k + l)y(k + l) 

+4>*(O k , k)H*(k)m - U d (k)r](k) (3.12) 


Thus with W(k) = n m (k) + 6, and using Eq. (3.12), the modal spatial velocity V m (k ) G 'FlFW 
for the k ih body is given by 


V m (k) ^ 


( 


ri(k) 

V{k) ) 


$*{k + 1, k)V m (k + 1) + n*(k)x(k) € 


(3.13) 


where the interbody transformation operator $(., .) and the modal joint map matrix 'H(k) are 
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defined as 


*(* + 1,*) = 


H(k) ± 


( 


o [n l (k + i)}*<f>{t k+ 1 ,k) 


\ 


^J7(k+i)x77(k) 


0 <j>(k + 1, k) J 

\ 


i -[n d (fc)]- 

^ 0 Hyr(k) ) 


n X(k)xJ?(k) 


where 


Note that 


where 


■*(*) = 


Hr{k)t H{k) 4 >{O k ,k)en n rW 


X6 


$(k + l,k) = A(k + l)B{k + l,k) 


/ 


[n *(*)]• 

<f>(k,tk) 


e n W(k)x6 and B{k + i,k) = [o, <f>{t k+1 , k)} £ n 6xWw 


The modal joint map matrix 'H(k) can be partitioned as 


H(k) = 


' «/(*) N 


n r (k) 


j^M(k)x77(k) 


where 


H f (k) = [/, — [n d (A;)]“] € and H r {k) = [0, H{k)<f>(O k ,k)\ £ K n ’ 


(3.14) 

(3.15) 

(3.16) 

(3.17) 

(3.18) 

(fc)xA^(fc) 

(3.19) 
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With W — J2kL x AT(Ar), define the spatial operator £$ as 




0 

*(2,1) 

0 


0 

0 

$(3,2) 


0 

0 

0 


0 

0 

0 


0 ... $(N,N- 1) 0 ) 


€ 71 ' 


A Tx77 


(3.20) 


Noting that £<$ is nilpotent (i.e. £$ = 0), we define the spatial operator $ as 


( 

$ = [/ - f*]" 1 = / + £* + ••• + = 

V 

(3.21) 

where 

$(i,j) = $(v - 1) • • • $(j + 1 ,j) for i > j 

Also define the spatial operator Ti = diag^7-^( A:) E 7 . Using these spatial operators, 

and defining V m = col|V m (fc)} E 7Z^, from Eq. (3.13) it follows that the spatial operator 
expression for V m is given by 


I 

$( 2 , 1 ) 


0 

/ 


0 

0 


e n 


77x77 


$(7V, 1) $(AT,2) ... I J 


V m = $*ft*X (3.22) 

3.2 Modal Mass Matrix for a Single Body 

In this section we derive an expression for the modal mass matrix of the k th body. With 
Vs(jk) € 'R 6 denoting the spatial velocity of node jk, and V s (k) = col|V' s (j* ; ) j- 6 the 

vector of all nodal spatial velocities for the k th body, it follows (as in Eq. (3.7)) that 

V 3 (k) = B*(k)V(k) + u(k ) = [11(A), B*(k)]V m (k) (3.23) 
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where 


B(k) = [^(fc,U),^(A;,20,-'-,^,H s (fc))]G^ 6x6 "^ 


(3.24) 


Since M s (k ) is the structural mass matrix of the k th body, the kinetic energy of the k th body 
is given by 

\v: (k)M,(k)V.(k) = \vi(k)M m (k)V m (k) 


where 





A 

' n*(fc) ' 

1 .... . 

. 1 

' n*(k)M s {k)n{k) n *(k)M s (k)B*(k) ^ 

M m {k) = 


B\k ) j = 



k B(k) , 

I 

! 

v B{k)M,{k)U{k) B(k)M a (k)B*(k) j 


' Mink) 

Mint) ^ 

£ ^(fc)xAT(fc) 

(3.25) 


K M’Jik) 

MX(k) j 




Corresponding to the generalized speeds vector x(^)? M m (k) as defined above is the modal 
mass matrix of the k ih body. In the block partitioning in Eq. (3.25), the superscripts / 
and r denote the flexible and rigid blocks respectively. Thus M£f(k) represents the flex/flex 
coupling block, while Mj£(k) the flex/rigid coupling block of M m {k). We will use this 
notational convention all through this publication. Note that M™(k) is precisely the rigid 
body spatial inertia of the k th body. Indeed, M m (k) reduces to the rigid body spatial inertia 
when the body flexibility is ignored, he., no modes are used, since in this case n m (k) = 0 
(and so TI( A:) is null). 

Since the vector l(k,jk) from Tk to node depends on the deformation of the node, 
the operator B(k) is also deformation dependent. From Eq. (3.25) it follows that unlike the 
block M{1 (A), which is deformation independent, both the blocks M£(k) and (k) are 
deformation dependent. The detailed expression for the modal mass matrix can be defined 
using modal integrals which are computed as a part of the finite-element structural analysis 
of the flexible bodies. Such an expression for the modal mass matrix of the k th body is given 
in Appendix B in Eq. (B.8). Often the deformation dependent parts of the modal mass 
matrix are ignored, and free-free eigen-modes are used for the assumed modes TI ( A: ) . When 
this is the case, M^{k) is zero and M"(fc) is block diagonal. 
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3.3 Recursive Propagation of Accelerations 


Differentiating the velocity recursion equation, Eq. (3.13), we obtain the following recursive 
expression for the modal spatial acceleration ot m (k ) £ 7Z ^ for the k th body: 


a m (k) = V m (k) = 


^ v( k ) ^ 


cx{k) 


$*(& + 1, k)a m (k + 1) + W*(A:)x(fc) + a m{k ) 


( 3 . 26 ) 


where a(k) = V(k), and the Coriolis and centrifugal acceleration term a m ( k ) £ is 

given by 


a m (k) 


d$*(& + 1) k) 

dt 


V m {k + 1) + 


dtt»(fc) 

dt 


X( k ) 


( 3 . 27 ) 


The detailed expressions for a m (k ) can be found in Eq. (B.ll) in Appendix B. Defining 
a m = co l{«m(^)} £ and a m = col|a m (A:)| £ BlF , Eq. (3.26) can be reexpressed using 
spatial operators in the form 


o m = $*(7f*X + a m) 


( 3 . 28 ) 


The vector of spatial accelerations of all the nodes for the k th body, a 3 (k ) 
£ 7£ 6n »(*) } is obtained by differentiating Eq. (3.23): 


col{a s (i*:)} 


where 


a,( k) = V,(k) = [n (k), B*(k)]a m { k) + a(k) 


( 3 . 29 ) 


a(Ar) ^ col{a(j*)} = ^^^ V m (k) £ 7l 6n ^ 


( 3 . 30 ) 


3,4 Recursive Propagation of Forces 

We now develop the equations of motion for the k ih body. Let f(k — 1) £ 1Z 6 denote the 
effective spatial force of interaction, referred to frame Tk -\ , between the k th and 
( k — 1)^ bodies across the (k — l) ih hinge. Recall that the ( k — l) f/l hinge is between node t * 
on the k th body and node d^-i on the (k — l) th body. With f s (jk) £ 7£ 6 denoting the spatial 
force at a node j k , the force balance equation for node t k is given by 

f s (tk ) = 4>{tk,k - l)f(k - 1 ) + M s (t k )a s (t k ) + b(t k ) + /*(«*) ( 3 . 31 ) 
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For all nodes other than node t k on the k th body, the force balance equation is of the form 


fs(jk) - M a (jk)oc a (j k ) + b(j k ) + 


(3.32) 


In the above, fx{k) = K a (k)u(k) G 1l 6n d k ) denotes the vector of spatial elastic strain forces 
for the nodes on the k th body, while b(j k ) G K 6 denotes the spatial gyroscopic force for node 
jk and is given by 


KJk) = 


( 


u{jk)J{jk)u(jk) 

m{jk)u(jk)u{jk)p{jk) 


\ 


e it 


(3.33) 


where u)(j k ) € TZ 3 denotes the angular velocity of node j k . Collecting together the above 
equations and defining 


C(k,k-1) = 


<t>(t k ,k- 1) 


€ and b{k) = col [b{j k )} G 7^ 6n ‘ w (3.34) 


it follows from Eq. (3.31) and Eq. (3.32) that 

f a (k) = C(k, k - l)f(k - 1) + M s (k)ce s (k) + b(k) + K a (k)u(k) (3.35) 


where f s (k) = col|/ s (jt)| G 1Z 6n, ^ k \ Noting that 

m = B(k)f a (k) 


(3.36) 


and using the principle of virtual work, it follows from Eq. (3.23) that the modal spatial 
forces fm(k) G for the k ih body are given by 


fm(k) = 


/ n*(fc) ^ 


m 


/.(*) = 


‘ n •(*)/.(*) x 


V 


m 


(3.37) 


1 
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Premultiplying Eq. (3.35) by 


and using Eq. (3.25), Eq. (3.29), and Eq. (3.37) 


1 n-(ir) ' 


B(k) 


leads to the following recursive relationship for the modal spatial forces: 


/»(*) = 


( n*{k)C(k,k-i ) N 


f(k - 1 ) + M m (k)a m {k ) + b m {k ) + I< m (k)d(k) 


1 [n‘(*r))- N 


^ B(k)C(k,k- 1 ) j 

4>{t k , k - 1 )f(k - 1 ) + M m (k)a m (k) + b m (k ) + K m (k)d(k) 

( 3 . 38 ) 


\ <i>{k,t k ) J 

= <!>(&, k - 1 )fm(k - 1 ) + M m (k)a m (k) + b m (k) + I< m (k)$(k) 


Here we have defined 


b m {k) 


and the modal stiffness matrix 


K m (k) = 


n*(t) ' 

B(k) j 


[ 6 (Jt) + M,(k)a(k)) e 


( 3 . 39 ) 


Tl*(k)K a (k)U(k) 0 
0 0 




( 3 . 40 ) 


The expression for K m (k) in Eq. (3.40) uses the fact that the columns of B*(k ) are indeed the 
deformation dependent rigid body modes for the k th body and hence they do not contribute 
to its elastic strain energy. Indeed, when a deformation dependent structural stiffness matrix 
K s (k) is used, we have that 


K 3 (k)B*(k) = 0 ( 3 . 41 ) 

However the common practice (and also followed here) of using a constant, deformation- 
independent structural stiffness matrix leads to the anomalous situation wherein Eq. (3.41) 
does not hold exactly. In view of this, without any loss in accuracy we ignore these anomalous 
extra terms and drop them from the left-hand side of Eq. (3.41). 

The velocity-dependent bias term b m (k) is formed using modal integrals generated 
by standard finite-element programs, and a detailed expression for it is given in Eq. (B.50) 
in Appendix B. From Eq. (3.38), the operator expression for the modal spatial forces 
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(3.42) 


f m = col|/ m (fc)| G lF for all the bodies in the chain is given by 

f m — + b m + A m $) 

where 

M m = diag{M m (fc)} € iF * 17 , I< m = diag{/F m (*)} € and b m = col{6 m (fc)} € lF 

From the principle of virtual work, the generalized forces vector T G RF for the multibody 
system is given by the expression 

T = HU (3.43) 


3.5 Operator Expression for the System Mass Matrix 

Collecting together the operator expressions in Eq. (3.22), Eq. (3.28), Eq. (3.42) and Eq. (3.43), 
we obtain the following: 

V m = H* X 

a m = $*(«** + ««) (3.44) 

fm - HM m a m + b m + K m d) = $C**«*x + $(Mm$*a m + b m + I< m d) 

T = Hf m = H<S>M m <f>-H*x + H$(M m $*a m + b m ) 

= Mx + C 

where 

M = H$M m $*H* G and C = H$(Mm$*a m + b m + K m $) G RF (3.45) 

Here A4 is the system mass matrix for the serial chain and the expression is 

referred to as the Newton-Euler Operator Factorization of the mass matrix. The term C is 
the vector of Coriolis, centrifugal, and elastic forces for the system. 

A noteworthy fact about the operator expressions for M and C is that they are iden- 
tical in form to the corresponding expressions for rigid multibody systems (see references 
[1, 5]). Indeed, the similarity is more than superficial, and the key properties of the spatial 
operators that are used in the analysis and algorithm development for rigid multibody sys- 
tems also hold for the spatial operators defined here for flexible multibody systems. Apart 
from the pedagogical importance, a significant advantage of this is that a large part of the 
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analysis and algorithms for rigid multibody systems can be easily carried over and applied 
to flexible multibody systems. This is precisely the approach we adopt here. 


4. Inverse Dynamics Algorithm 


This section describes a recursive Newton-Euler inverse dynamics algorithm for computing 
the generalized forces T, for a given set of generalized accelerations x an d system state 
The inverse dynamics algorithm also forms a part of forward dynamics algorithms 
such as those based upon composite body inertias or the conjugate gradient method [6]. 

Collecting together the recursive equations in Eq. (3.13), Eq. (3.26), Eq. (3.38) and 
Eq. (3.43) we obtain the following recursive Newton-Euler inverse dynamics algorithm: 


V m (iV + l) = 0, <x m (N + 1) = 0 

for k = N • ■ ■ 1 

< V m (k) = $*(k + l,k)V m (k + l) + H-(k) x (k) 

a m (k) = $*(k+l,k)a m (k + l) + H*(k)x(k)+a m (k) 

end loop 

( 4 . 1 ) 


fm( 0) = 0 

for k = 1 • • • N 

< fm(k) - $(k, k - 1 )f m (k - 1) + M m (k)a m (k) + b m (k) + K m (k)ti(k) 

T(k ) = n(k)f m (k) 

end loop 

The structure of this algorithm closely resembles the recursive Newton-Euler inverse dynam- 
ics algorithm for rigid multibody systems [7, 1]. It assumes without any loss in generality 
that the base body is stationary and also that the tip force is zero. Base mobility is handled 
by attaching an additional 6-dof hinge between the mobile base and an inertial frame. In 
case there is a non-zero tip force, / m ( 0) should be initialized to the value of the tip force 
in the algorithm. By taking advantage of the special structure of <&(k + l,k) and 'H(k) in 
Eq. (3.14) and Eq. (3.15), the Newton-Euler recursions in Eq. (4.1) can be further simplified. 
Using the superscripts / and r as before to denote the flexible and rigid components, we 
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have the following partitioning: 


vm ) ... 

' °L(k) ^ 

. Uk) = 

' fm(k) 

, and T(k) = 

' T'(k ) 'l 

KTO ) 

{ <X r m( k ) ) 


\ fm t 


V T’(k) j 


Based upon this partitioning, the simplified inverse dynamics algorithm is as follows: 


V m (N + 1) = 0, a m (N + 1) = 0 
for k = N • • • 1 

vm = m 

Vm(k) ~ <l>*(tk+uk)’A*(k + l)Kn(A: + 1) + Hjr(k)[3(k) — U d (k)ri(k) 

a L( k ) = TO 

a m(k) = <^*(^Jk+i5 k)A*(k + l)a m (fc + 1) T H^r{k)^(k) — + a m( k ) 

end loop 

( 4 . 2 ) 


U 0) = o 

k = 1 -N 

1 )fmi k ~ 0 d" + b m (k) + K m (k)d(k) 

' sm - [imm« ^ 

k HAk)Uk) ) 

Flexible multibody systems have actuators typically only at the hinges. Thus for the 
k th body, only the subset of the generalized forces vector T(k ) corresponding to the hinge 
actuator forces T r (k) can be set, while the remaining generalized forces T*(k) are zero. Thus 
in contrast with rigid multibody systems, flexible multibody systems are under- actuated sys- 
tems [8], since the number of available actuators is less than the number of motion dofs in 
the system. For such under- actuated systems, the inverse dynamics computations for the 
generalized force T are meaningful only when the prescribed generalized accelerations \ form 
a consistent data set. For a consistent set of generalized accelerations, the inverse dynamics 


/m(fc) = A(k)<f>(tk,k - 
T(k ) = 


1 T f (k) ' 


\ T r (k) 


end loop 
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computations will lead to a generalized force vector T such that T*(.) = 0. 


5. Composite Body Forward Dynamics Algorithm 


The forward dynamics problem for a multibody system consists of computing the generalized 
accelerations x for a given set of generalized forces T and state of the system The 

composite body forward dynamics algorithm consists of the following steps: (a) computing 
the system mass matrix M , (b) computing the bias vector C, and (c) numerically solving 
the following linear matrix equation for \\ 

Mx — T- C (5.1) 

Later in Section 6 we describe the recursive articulated body forward dynamics algorithm 
that does not require the explicit computation of either A4 or C. 

It is evident from Eq. (5.1) that the components of the vector C are the generalized 
forces for the system when the generalized accelerations x are a U zero. Thus C can be 
computed using the inverse dynamics algorithm in Eq. (4.2). An efficient composite-body- 
based recursive algorithm for the computation of the mass matrix M is described in this 
section. This algorithm is based upon the following lemma which describes a decomposition 
of the mass matrix into block diagonal, block upper triangular and block lower triangular 
components. 


Lemma 5.1: Define the composite body inertias R(k) E R?* 7(k)xJT(k) recurs i ve ly for all 

the bodies in the serial chain as follows: 


R( 0) = 0 
for k = 1 • • • N 

(5.2) 

R(k) = $(k,k-\)R{k-l)$*{k,k-l) + M m {k) 

end loop 


Also define R = diag{-R(fc)| £ . Then we have the following spatial operator decom- 
position where $ = — I: 


$M m $* = R + $R + R$* 


(5.3) 
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Proof: See Appendix A. 


Physically, R(k) is the modal mass matrix of the composite body formed from all 
the bodies outboard of the k th hinge by freezing all their (deformation plus hinge) dofs. It 
follows from Eq. (3.45) and Lemma 5.1 that 

M = = HRH* + HQKHT + HR&H* ( 5 . 4 ) 

Note that the three terms on the right of Eq. (5.4) are block diagonal, block lower triangular 
and block upper triangular respectively. The algorithm for computing the mass matrix 
M computes these terms recursively. The main recursion proceeds from tip to base, and 
computes the blocks along the diagonal of A4. As each such diagonal element is computed, 
a new recursion to compute the off-diagonal elements is spawned. Its structure is similar 
to that of the composite body algorithm for computing the mass matrix of rigid multibody 
systems [6, 9], and is as follows: 


R{ 0 ) = 0 

for k = 1 • • • N 

R(k) = <f>(k,k-l)R(k-l)$*{k,k-l) + M m {k ) 

= A(k)<f>(t k , k - 1 )R rr (k - 1 )<j>*{t k , k - I )A'{k) + M m {k) 

X(k) = R{k)H*{k) 

M s (k, k) = H{k)X{k) 

< , 

for j = (k + 1) • • • N 

XU) = m,i - i)x(j - 1 ) = - 1 ) 

M(j,k) = M*(k,j) = 'HU)X(j) 

end loop 

end loop 

( 5 . 5 ) 

The structure of the above algorithm for computing the mass matrix closely resembles the 
composite rigid body algorithm for computing the mass matrix of rigid multibody systems 
[6, 9]. Using the sparsity of both 'Hf(k) and 'H T {k), additional computational simplification 
of the steps in the above algorithm is easy to implement. 
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6. Factorization and Inversion of the Mass Matrix 


An operator factorization of the system mass matrix M, denoted the Innovations Operator 
Factorization, is derived in this section. This factorization is an alternative to the Newton- 
Euler factorization in Eq. (3.45) and, in contrast with the latter, the factors in the Innovations 
factorization are square and invertible. Operator expressions for the inverse of these factors 
are developed and these immediately lead to an operator expression for the inverse of the 
mass matrix. Using further operator identities, we obtain an operator expression for the 
generalized accelerations x in terms of the applied generalized forces T . The recursive 
implementation of this expression leads to the recursive articulated body forward dynamics 
algorithm described in Section 7. The operator factorization and inversion results here 
closely resemble the corresponding results for rigid multibody systems (see [1]). 

Given below is a recursive algorithm which defines some required articulated body 
quantities: 


< 


P*( 0) = 0 


• k = 1-N 


P(k) = 

$(] k,k - l)P + {k- l)$*(Jb,lfc- 1) 

D(k) = 

H{k)P{k)H*{k) G 

cm = 

P(k)H*(k)D-\k) G 

K(k + 1, k) — 

$(* + 1, *)<?(*) G 

f(k) = 

I - G{k)H{k) G 

P+(k) = 

~(k)P(k) G 

+ 1, k) = 

$(k + l,k)T{k) G 7^(*)*V(fc) 


+ M m (k) G 7 ^U)xV(fc) 


end loop 


( 6 . 1 ) 


The operator P G } s defined as the block diagonal matrix with the k th diagonal 

element being P(k). The quantities defined in Eq. (6.1) form the component elements of the 
following spatial operators: 

D = TiPH* = diag{£>(ifc)} 6 K*** 

G = PWD = diag{G(fc)} € 1l Wx " 

I< = 6<tG G n WxAf 
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r = I -GH = diag{r(fc)} e 

£* = EtTC'R** 77 ( 6 . 2 ) 

The only nonzero block elements of K and Eq are the elements’ Jlf(fc + 1, k)'s and *(& + 1, fc)’s 
respectively along the first sub-diagonal. 

As in the case for £$, Ey is nilpotent, so we can define 

I 0 ... 0 

*(2,1) I ... 0 

*(AT,1) *(iV, 2) ... I 

where 

*(i,j) = *(i,i - 1) *(i + l,j) for i > j 

The structure of the operators Ey and * is identical to that of the operators E$ and * 
respectively except that the component elements are now *(i,j) rather than *( 1 , 7 ). Also, 
the elements of * have the same semigroup properties as the elements of the operator 3>, 
and as a consequence, high-level operator expressions involving them can be directly mapped 
into recursive algorithms, and the explicit computation of the elements of the operator * is 
not required. 

The Innovations Operator Factorization of the mass matrix is defined in the following lemma. 

Lemma 6.1: 

M = [I + H$I<]D[I + H$I<}* (6.4) 

Proof: See Appendix A. 1 

Note that the factor [I + 7{$K] £ 7 ^*^ j s square, block lower triangular and non- 
singular, while D is a block diagonal matrix. This factorization may be regarded as a block 
LDL* decomposition of M. The following lemma gives the closed form operator expression 
for the inverse of the factor [I + H<bK]. 
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Lemma 6.2: 


[/ + mK]~ l = [i - mi<] ( 6 . 5 ) 

Proof: See Appendix A. ■ 

It follows from Lemmas 6.1 and 6.2 that the operator expression for the inverse of the mass 
matrix is given by: 

Lemma 6.3: 

AT 1 = [/ - WlK]*D- l \I - H9K] (6.6) 


Once again, note that the factor [I — H9K] is square, block lower triangular and 
nonsingular and so Lemma 6.3 may be regarded as providing a block LDL * decomposition 
of M~\ 


7. Articulated Body Forward Dynamics Algorithm 


Wc first use the operator expression for the mass matrix inverse developed in Section 6 to 
obtain an operator expression for the generalized accelerations x • This expression directly 
lea,ds to a recursive algorithm for the forward dynamics of the system. The structure of 
this algorithm is completely identical in form to the articulated body algorithm for serial 
rigid multibody systems. The computational cost of this algorithm is further reduced by 
separately processing the flexible and hinge dofs at each step in the recursion, and this forms 
the articulated body forward dynamics algorithm for serial flexible multibody systems. 

The following lemma describes the operator expression for the generalized accelera- 
tions x * n terms of the generalized forces T. 


Lemma 7.1: 

X = [/ - HVKYD-' [T - HV{KT + Pa m + b m + I< m d}] - K*V*a m (7.1) 
Proof: See Appendix A. ® 
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As in the case of rigid multibody systems [1, 3], the direct recursive implementation 
of Eq. (7.1) leads to the following recursive forward dynamics algorithm: 


z+( 0) = 0 


for k = 1 • • • n 


z(k) 

e(k) 

v{k) 

z+ (k) 


$(k, k - \)z+(k - 1) + P(k)a m (k) + b m (k ) -f I< m (k)'d(k ) 

T(k) - H(k)z{k) ( 7 - 2 ) 

D~ 1 (k)e(k) 

z{k) + G(k)e{k) 


[ end loop 


&m{ n + 1 ) ~ 0 

for k = n •• • 1 

< 4 (*) = + 1 ) 

X(k) = v(k) - G'(k)a*(k) 
a m (k) = a t,(k) + H-(k)x(k) + a m (k) 

end loop 

All the dofs for each body as characterized by its joint map matrix H*{.) are processed to- 
gether at each recursion step in this algorithm. However, by taking advantage of the sparsity 
and special structure of the joint map matrix, additional reduction in computational cost is 
obtained by processing the flexible dofs and the hinge dofs separately. These simplifications 
are described in the following sections. 


7.1 Simplified Algorithm for the Articulated Body Quantities 

We describe intuitively the basis for the separation of the modal and hinge dofs for the body. 
First we recall the velocity recursion equation in Eq. (3.13) 

V m (k) = + 1 ,k)V m {k + 1) + H*(k)x(k) (7.4) 
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and the partitioned form of h((k) in Eq. (3.15) 


n(k) = 


' H ,(k) ^ 

K Hr (*) j 


Introducing a dummy variable k', we can rewrite Eq. (7.4) as 


v m (k') = <S>*(k + l,k')V m (k+l) + H}(k)ri(k) 
V m (k) = $*(k',k)v m (k') + n;(k)/3(k) 


( 7 . 5 ) 


( 7 . 6 ) 


where 


$(fc + 1, k') = ${k + 1, k ) and $(fc', k) = I 

Conceptually, each flexible body is now associated with two bodies. The first one has the 
same kinematical and mass/inertia properties as the real body and is associated with the 
flexible dofs. The second body is a fictitious body and it is massless and has zero extent. It 
is associated with the hinge dofs. The serial chain now contains twice the number of bodies 
as the original one, with half the new bodies being fictitious ones. The new 7 -f* operator will 
have the same number of columns but twice the number of rows as the original 7i* operator. 
The new $ operator will have twice the number of rows as well as columns as the original one. 
Going through the same analysis as described in the previous sections, we once again obtain 
the same operator expression as Eq. (7.1). This expression also leads to a recursive forward 
dynamics algorithm as in Eq. (7.3). However each sweep in the algorithm will contain twice 
as many steps as the original algorithm. However each step will be processing only a smaller 
number of dofs leading to a reduction in the overall cost. Note that a new set of articulated 
body quantities also needs to be defined using an algorithm whose structure is similar to 
that of the algorithm in Eq. (6.1). 

For low-spin multibody systems, an additional source of computational simplification 
is possible. Given the inherent linearization that results from using modes for modeling body 
deformation, there is typically little loss in model fidelity for such systems when the defor- 
mation and deformation rate dependent terms in M and C are dropped from the dynamical 
equations of motion [10], Such models have been dubbed ruthlessly linearized models. The 
ruthlessly linearized models are considerably less complex, and their use results in a sub- 
stantial reduction in the computational complexity of the dynamical algorithms. For these 
models, the modal mass matrices for the component bodies are constant in the body frame. 
In this case, given matrices A, B and C, the following matrix identity 

[A + BCB*]~ X = A~ x - A~ X B[C~ X + B* A~ x B)~ x B* A~ x ( 7 . 7 ) 
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can be used to simplify the formation of D~ l in Eq. (6.1). As a result of this simplification, 
the computational complexity of each recursion step for the articulated body inertias becomes 
a quadratic rather than a cubic function of the number of modes. 

Separating the flexible and hinge dofs for each body, and using the simplifications 
from Eq. (7.7) for a ruthlessly linearized model, a simplified recursive algorithm for the com- 
putation of the articulated body quantities is described below. When a ruthlessly linearized 
model is not desired the only change required in the algorithm below is that instead of com- 
puting indirectly, it must be computed by directly inverting the matrix Dj(k). The 

following matrices can be used only for the ruthless model and need to be precomputed just 
once prior to the dynamical simulation: 


for k = 1 • • • N 

m = \H } {k)M m {k)H){k)}-' en"*" 

<( k ) = / Hj{k)A{k) € TZ^ xG ( 7 - 8 ) 

T(Jfc) = A (k)((k) e K** 6 

n(k) = <*(*)T(fc) € n 6 * 6 

end loop 
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The remainder of the algorithm for computing the articulated body quantities is as follows: 


P+(0) = 0 
for k = 1 • • • N 

T(k) = B(k,k-l)P + {k-l)B*(k,k-l)eH 6x6 
p(k) = A(k)r(k)A*(k) + M m (k)enF^ x ^ 

D f (k) = Hj{k)P{k)H){k) € ^ nm(fc)xn '" (fc) 

Dj\k) = a (k) - T(it)[r- 1 (^) + n(k)]- l (k)r*(k) e 7e nm(fc)xnm(fc) 

G/(lfc) = P{k)n}{k)D]\k)eK 77 W* n ”'W 
rj(k) = / - G,(k)H,(k) e 

P r (fc) = T}{k)P{k) G 

D r (jfe) = K(fc)P r (ik)W:(fc) € ^ nr(fc)x " r(fc) 

G>(fc) = P r (ifc)^:(^)^ r _1 (^) e ^ Wxnr(fc) 

r r (Ar) = I -G r {k)H T (k) 

P+(k ) = T T (k)P r (k) G 

^(k + l,fc) = $(Jfc + 1, fc)r(fc) G 7 

end loop 

The sparsity of B(k-\- 1, jfc), 'Hf(k) and H T (k) leads to still further simplification of the above 
algorithm. Using the symbol “x” to indicate “don’t care” blocks, the structure in block 
partitioned form of some of the quantities in Eq. (7.9) is given below: 

T (k) = <f>(tk,k - l)Pfi(k — l)<f>*(tk,k — l), (P^(fc) is defined below) 

x \ 

Gf(k) — , where g(k) = fi(k)Dj 1 (k) G 7£ 6xnm U) 

y 9(k) j 

and fi(k) = [P rf (k), P TT (k)]H m j{k ) G 7l 6xn ">^ 
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p r (k) = 


X X 


\ 


X P R (k) 


, where P R (k) = P rr {k ) - g(k)f(k) E 71' 


6X6 


D r (A:) = Hr{k)P R {k)H}(k)<=1l n ' ik)xnrW 


G r (k) = 


T r (^) = 


P + (i) = 


/ \ 

X 


Gr(*) 


/ x 

0 Tr(£) 


, where <?*(&) = Pr^T/X^LT 1 ^) € 7l 6xn '<*> 

\ 


, where r R (k) = 1 — G R {k)Hjr(k ) E 7?. 


6X6 


'* X ' 


V 


, where Pr(^) = r R {k)P R {k ) E 7£ 


6x6 


y 


X Pr^) 

Using the structure described above, the final simplified algorithm for computing the artic- 
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ulated body quantities is as follows: 


p+(0) = 0 

for k = 1 • • • N 

r (k) = <t>(t k ,k-i)p£(k-i)p(t k ,k-i) 

P(k ) = A(k)T(k)A*(k) + M m (k) 

Dj(k) = Hj(k)P(k)H}(k) 

D]\k ) = A(ifc) — TCifcjtr-'CJb) + n(jfc)]- 1 (fc)'r*cifc) 

p(k) = [ P rf (k ), p~(*)]w;(*) 

g(k) = n{k)Dl l {k) 
p R (k) = P rr (k) - <?(%*(&) 

Gfl(fc) = PR(k)H*r(k)Dj(k) 
r R (k) = I-G R {k)Hr{k ) 

P+(k) = T R (k)P R (k) 

end loop 


( 7 . 10 ) 


7.2 Simplified Articulated Body Forward Dynamics Algorithm 


The simplified recursive articulated body forward dynamics algorithm for a serial flexible 
multibody system follows directly from the recursive implementation of the expression in 
Eq. (7.1) for the multibody system with fictitious bodies described in Section 7.1. It consists 
of the following steps: (a) computation of the articulated body quantities using Eq. (7.8) 
and Eq. (7.10), (b) a base-to-tip recursion as in Eq. (4.2) for computing the modal spatial 
velocities V m (k), and the bias terms a m (k ) and b m (k ) for all the bodies, and (c) a tip-to-base 
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recursion followed by a base-to-tip recursion for the joint accelerations x as described below: 


4(0) — o 


for k = 1 • • • N 

*(*) = 


f (lr\ ^ 

Zf{k) 


Zr(k) 


= A{k)<t>(t k , k - 1 )4(* - 1) + 


ps r (k) 

P rr {k) 

+b m (k) + K m (k)ti(k) e 


\ 


a m R{k :) 


/ 


e f {k) = 

T f (k) - z f (k) + [U d (k)]*z r (k) 

v s {k) = 

Dj\k)e } (k) E n nm W 

z R (k) = 

z T (k) + g(k)e f (k)en 6 

efl(*0 = 

T R (k) - Hf(k)z R (k) E 

v R (k) = 

D R l {k)e R {k) E 1l n ^ 

4(*0 - 

z R {k) + G R (k)e R (k) E K 6 

end loop 



( 7 . 11 ) 


a m(^ + 1) — 0 

for k = N • • • 1 

a£(fc) = 4>*(t k+u k)A m (k + l)a m (k + 1 ) e n 6 
m = u R (k) - GUk)aUk) e ft nr W 

a R (k) = a^(k) + H^(k)0{k) + a mR {k)eH 6 

ij(k) = u } (k) - g‘(k)a R (k) E 

( \ 

G 7 Z&W 


a m (k) = 

end loop 


ij(k) 

a R (k) - n d (k)rj(k) ) 


( 7 . 12 ) 


Simplifications arising from the sparsity of 'H(k), along the lines described in Section 7.1, 
have been incorporated in the above algorithm. In contrast with the composite body forward 


33 



dynamics algorithm described in Section 5, this algorithm does not require the explicit 
computation of either M, or C. The structure of this articulated body algorithm closely 
resembles the recursive articulated body forward dynamics algorithm for rigid multibody 
systems described in references [11, 1]. 


8. Computational Cost 


This section discusses the computational complexity of the composite body and the articu- 
lated body forward dynamics algorithms* 

Flexible multibody systems typically involve both rigid and flexible bodies and, in 
addition, different sets of modes are used to model the flexibility of each body. As a con- 
sequence, where possible, we describe the contribution of a typical (non-extremal) flexible 
body, denoted the k th body, to the overall computational cost. Note that the computational 
cost for extremal bodies as well as for rigid bodies will typically be lower than that for a 
non-extremal flexible body. Summing up this cost for all the bodies in the system will give 
a figure close to the true computational cost for the algorithm. Without any loss in gener- 
ality, we have assumed here that all the hinges are single dof rotary joints and that free-free 
assumed modes are being used. 

All costs given below are based on the use of the ruthlessly linearized dynamics model 
of the flexible body [10, 12] wherein all deformation and deformation rate dependent terms 
are dropped from Ai and C. It has been pointed out in recent literature [13, 10] that the use 
of modes for modeling body flexibility leads to “premature linearization” of the dynamics, 
in the sense that while the dynamics model will contain deformation dependent terms, the 
geometric stiffening terms will be missing. Indeed, these missing geometric stiffening terms 
are typically the dominant first-order (deformation) dependent terms. When the body 
spin rates are high, it is necessary to take additional steps to include the stiffness terms 
to obtain a “consistently” linearized model with the proper degree of fidelity. However for 
systems with low spin rate, the contribution of the first order terms to the dynamics model is 
negligible. Dropping the deformation 7](k) and deformation rate rj dependent terms leads to 
the “ruthlessly” linearized dynamics model and significant reductions in the computational 
complexity of the algorithms. In the ruthlessly linearized model, Eq. (B.8), Eq. (B.ll) and 
Eq. (B.50) are approximated as follows: 


M m (k) » M°(Ar), a m (k) 


l mR 


(k) 


and b m (k ) ss b° m (k) 


( 8 . 1 ) 


Note that with this approximation, M m (k ) is constant in the body frame, while a m (k) and 
b m (k) are independent of Tj(k) and fj(k). A large number of flexible multibody systems 
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are indeed low spin systems for which the ruthlessly linearized model is adequate. The 
computational costs described below are for such a system and are based on the use of a 
ruthlessly linearized model. 


8.1 Computational Cost of the Composite Body Forward Dy- 
namics Algorithm 

The composite body forward dynamics algorithm described in Section 5 is based on solving 
the linear matrix equation 


Mx=T-C 


The computational cost of this forward dynamics algorithm is given below: 


1. Cost of computing R(k ) for the k th body using the algorithm in Eq. (5.5): 

products = 48n m (A;)+90 

97 

additions = n^(k) —n m (k) 116 

2. Contribution of the k th body to the cost of computing M (excluding cost of R(k)'s) 
using the algorithm in Eq. (5.5): 

products = k\\2n 2 m (k) + 34n m (fc) + 13] 
additions = &[lln^(A:) + 24n m (k) + 13] 

3. Setting the generalized accelerations \ = 0> the vector C can be obtained by using the 
inverse dynamics algorithm described in Eq. (4.2) for computing the generalized forces 
T. The contribution of the k th body to the computational cost for C(k ) is 

products = 2n 2 m {k) + 54 n m (k) + 206 
additions = 2n^(fc) + 50 n m (k) + 143 

4. The cost of computing T — C is 


products = 0 
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additions = A r 

5. The cost of solving the linear equation in Eq. (5.1) for the accelerations x I s 

products = -A f 3 + -N 2 — —Af 
6 2 3 

additions = -J\f 3 + Af 2 — 7 Af 
6 0 

The overall complexity of the composite body forward dynamics algorithm is 0(Af 3 ). 

8.2 Computational Cost of the Articulated Body Forward Dy- 
namics Algorithm 

The articulated body forward dynamics algorithm is based on the recursions described in 
Eq. (7.8), Eq. (7.10), Eq. (7.11) and Eq. (7.12). Since the computations in Eq. (7.8) can be 
carried out prior to the dynamics simulation, the cost of this recursion is not included in the 
cost of the overall forward dynamics algorithm described below: 

1. The algorithm for the computation of the articulated body quantities is given in 
Eq. (7.10). The step involving the computation of D~ i (k ) can be carried out ei- 
ther by an explicit inversion of D(k ) with 0(n^(£)) cost, or by the indirect procedure 
described in Eq. (7.10) with 0(n^(A:)) cost. The first method is more efficient than 
the second one for n m (k) < 7. 

• Cost of Eq. (7.10) for the k th body based on the explicit inversion of D(k ) (used 
when n m (k) < 7): 

5 25 764 

products = ~n 3 m (k) + — n 2 m (k ) + -jj-n m (fc) + 180 

additions = |n^(fc) + — n^( k ) + —z~ri m (k) + 164 
o z 0 


• Cost of Eq. (7.10) for the k th body based on the indirect computation of D l (k) 
(used when n m (k) > 8): 

products = 12n^(fc) + 255 n m (k) + 572 
additions = 13n^(fc) + 182n m (fc) + 445 
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Number of Bodies 

Number of Modes n m 

N 

5 

10 

15 

3 

5223m + 4195a 

11432m + 9677a 

18587m + 16287a 

6 

12591m + 10075a 

27611m + 23315a 

44846m + 39225a 

9 

19591m + 15955a 

43790m + 36953a 

71105m + 62163a 


Table 1: Computational cost for the articulated body for- 
ward dynamics algorithm 


2 . 


The cost for the tip-to-base recursion sweep in Eq. (7.11) for the k th body is: 

products = n 2 m {k) + 25 n m (k) + 49 
additions = n^(fc) + 24n m (fc) + 50 


3. The cost for the base-to-tip recursion sweep in Eq. (7.12) for the k th body is: 


products = 18n m (A;) 4- 52 

additions = I9n m (k) + 42 


The overall complexity of this algorithm is 0(Nn 2 m ), where n m is an upper bound on the 
number of modes per body in the system. Thus it is to be expected that the articulated 
body algorithm will be more efficient than the composite body algorithm as the number 
of bodies increases. Indeed a comparison of the computational costs of the two forward 
dynamics algorithms reveals that the articulated body algorithm is the more efficient even 
when a small number of assumed modes is used. 

Table 1 describes the computational cost of the articulated body forward dynamics 
algorithm for a different number of bodies N and number of modes n m for a serial chain 
system. In the table, m and a denote floating point multiplications and additions respectively. 
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9. Extensions to General Topology Flexible Multi- 
body Systems 


For a rigid multibody system, reference [5] describes the extensions to the dynamics for- 
mulation and algorithms that are required as the topology of the system goes from a serial 
chain topology, to a tree topology and finally to a closed-chain topology system. The key 
to this progression is that the operator description of the system dynamics does not change 
despite the increase in the topological complexity of the system. Indeed, as seen here, the 
operator description of the dynamics remains the same even when the multibody system 
contains flexible rather than rigid component bodies. Thus, using the approach in [5] for 
rigid multibody systems, the dynamics formulation and algorithms for flexible multibody 
systems with serial topology can be extended in a straightforward manner to systems with 
tree or closed-chain topology. Rased on these observations, extending the serial chain dy- 
namics algorithms described in this publication to tree topology flexible multibody systems 
requires the following steps: 

L For any outward sweep involving a base to tip(s) recursion, at each body, the outward 
recursion must be continued along each outgoing branch emanating from the current 
body. 

2. For an inward sweep involving a tip(s) to base recursion, at each body, the recursion 
must be continued inwards only after summing up contributions from each of the 
incoming branches at the body. 


A closed-chain topology flexible multibody system can be regarded as a tree topology 
system with additional closure constraints. As described in reference [5], the dynamics 
algorithm for closed-chain systems consists of recursions involving the dynamics of the tree 
topology system, and in addition the computation of the closure constraint forces. The 
computation of the constraint forces requires the effective inertia of the tree topology system 
reflected to the points of closure. The algorithm for flexible multibody systems for computing 
these inertias is identical in form to the recursive algorithm described in [5]. 


10. Conclusions 


This publication uses the spatial operator algebra methodology to develop a new dynamics 
formulation and spatially recursive algorithms for flexible multibody systems. A key feature 
of the formulation is that the operator description of the flexible system dynamics is identical 
in form to the corresponding operator description of the dynamics of rigid multibody systems. 
A significant advantage of this unifying approach is that it allows ideas and techniques for 
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rigid multibody systems to be easily applied to flexible multibody systems. The algorithms 
use standard finite-element and assumed modes models for the individual body deformation. 

The Newton-Euler Operator Factorization of the mass matrix, M = 
forms the basis for recursive algorithms such as those for the inverse dynamics, the com- 
putation of the mass matrix, and the composite body forward dynamics algorithm for the 
flexible multibody system. Subsequently, the articulated body forward dynamics algorithm 
is developed, which, in contrast to the composite body forward dynamics algorithm, does 
not require the explicit computation of the mass matrix. The key sequence of steps in- 
volved in the development of the articulated body forward dynamics algorithm is: (a) 

the development of an alternative Innovations Operator Factorization of the mass matrix, 
M — [I + T{$K]D[I + (b) formation of the inverse of the factor [I + Ti^K]~ l = 

[I — TC^K], and (c) the formation of the operator expression for the mass matrix inverse, 
M~ x = [I — K]* D~ l [I — WflK], While the major focus in this publication is on the 

dynamics of flexible multibody systems with serial topology, the extension of the algorithms 
developed here to tree and closed chain topology systems is straightforward. 

Based on the discussion in reference [3] for rigid multibody systems, forward dynamics 
algorithms such as the conjugate gradient and triangularization algorithms can be extended 
to flexible multibody systems in a straightforward manner from the operator description 
of the dynamics presented here. While the computational cost of the algorithms depends 
on factors such as the topology and the amount of flexibility in the multibody system, in 
general, it appears that in contrast with the rigid multibody case, the articulated body 
forward dynamics algorithm is the more efficient algorithm for flexible multibody systems 
containing even a small number of flexible bodies. The variety of algorithms described 
here permits a user to choose the algorithm which is optimal for the multibody system at 
hand. The availability of a number of algorithms is even more important for real-time 
applications, where implementation on parallel processors or custom computing hardware is 
often necessary to maximize speed. 
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Appendix A: Proofs of the Lemmas 


At the operator level, the proofs of the lemmas in this publication are completely analogous 
to those for rigid multibody systems [1, 3]. 


Proof of Lemma 5.1: Using operators, we can rewrite Eq. (5.2) in the form 


M m — R — £<s,R£$ 


(A.l) 
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From Eq. (3.21) it follows that $£* = 5*$ = $ - 7 = 4>. Multiplying Eq. (A.l) from the 
left and right by $ and respectively leads to 


= $7?$* - <t>£*R£*<S>* = ($ + I)R($ + /)* - = R + &R + R$* 


Proof of Lemma 6.1: It is easy to verify that tPt* = tP. As a consequence, the recursion 
for P(.) in Eq. (6.1) can be rewritten in the form 

M m = P -£*P£l = P -£*P£* = P- £«P£% + K DK * (A.2) 

Pre- and post-multiplying the above by $ and respectively then leads to 

$M m $* = P + + 9KDKW 


Hence, 


=» M = = H\P + 3>P + P$* + *K DICV\K 

= D + HQKD + DK*$*H* + HQKDICW = [7 + H*K}D\I + mi<\' 


Proof of Lemma 6.2: Using a standard matrix identity we have that 

[7 + n^K]- 1 = 7 - n^[I + KH$]~ l K (A. 3) 


Note that 


'I'" 1 = 7 - £ q = (7 - £*) + £*GH = fc" 1 + I<7i (A A) 

from which it follows that 


tf- 1 # = 7 + I<H$ 

Using this with Eq. (A. 3) it follows that 

[7 + = 7 - miv- 1 $]-'k = 7 - m k 


■ 
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Proof of Lemma 7.1: From Eq. (3.44) and Eq. (3.45), the expression for the generalized 
accelerations \ is given by 


X = M~\T-C ) 

= [I-HVKYD- x [I-mK][T-H$[M m <t>*a m + b m + K m i)]] (A.5) 

From Eq. (A. 4) we have that 

[I - HVK)H$ = = m (A. 6) 

Thus Eq. (A.5) can be written as 

X = [I- H9KYD- 1 [T - m[KT + + b m + I< m d]\ (A.7) 

From Eq. (A. 2) it follows that 

M m = P-£*P£; => = + (A. 8) 

and so Eq. (A.7) simplifies to 

x = [/ - HVKYD- 1 [t - m[I<T + Pa m + b m + I< m d ] - HP$*a m ] (A.9) 

From Eq. (A. 4) we have that 

[7 - nyKYD-'HPQ* = [I - TWKYICV = 7f*^*[^"* - = ICW (A. 10) 

Using this in Eq. (A.9) leads to the result. I 


Appendix B: Expressions for M m (k ), a m {k ) and b m (k) 


The modal spatial displacement influence vector IT- 7 (7:) for node ji* has the structure: 


iP(Jfc) = 


' V(*) 
v iK k ) 


g ^6xn m (fc) 


(B.l) 


The components of the vectors A *(k) £ 7^ 3xnm (*) and ^(k) £ K 3x "' n W are the modal slope 
displacement influence vector and the modal translational displacement influence vector re- 


42 



spectively for node j k . They define the contribution of the various modes to the slope (or 
differential change in orientation) and translational deformation for the j\ [ h node on the 
k th body. Define 

Ujfc) = a j (k)r,(k) e n 3 , f>v(jk) = 7 j (k)m e n 3 , and 6 t (j k ) = 7 i (%(*) € K 3 (B.2) 
Note that 

l( k Jk) = lo {jk) + h{j k ) 

where l 0 (j k ) denotes the undeformed vector from frame T k to node j k . Recall from Eq. (3.1) 
that 

, ( Jtik) ™(h)p{jk) 

M s {jk ) = 

^ ~™{jk)p(jk) m(j k )I 

B.l Modal Integrals for the Individual Bodies 

Defined below are a set of modal integrals for the k th body which simplify the computation 
of the modal mass matrix M m {k ) and the bias vector b m (k). These modal integrals can be 
computed as a part of the finite-element structural analysis of the individual bodies. 

n .(*) 

m i k ) = £ m (i*) 

j~ik 

n ,(k) 

Po = [l/m(fc)] £ m (jk)[p(jk) + ?o( k ijk)] € 'll 3 

i=*k 

».(*) 

P k li r ) = [l/™(*0] £ ™(jkh’r( k ) € 

i= i* 

n,(t) 

£‘(r) s y. ™(h)bm - PVMm e ft 3 

j=ijt 

A _ . . 

F 0 k (r) = £ J(jk)Ki k ) + Jfc) + p(ife)]7r( fc ) - ™{jk)lo{ k ,jk)p(jk)K( k ) € K 

3 = Ifc 
"»(*) 

F\{ r i s ) = £ "*(iit)7i(^)[7r(*) - P(j*)^(*)] £ ft 3 

j=lfc 
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*»•(*) 

G k (r,s) = E [K( k )]*J(n)K( k ) + m Uk)lK( k )]*p(jkhl( k ) + m{j k )[\ 3 ,(k)Yp(j k )'y}(k) + 

:=ik 

m (jk)hr( k )]*li( k ) g TV 

Jo = - E JUk) - m(jk)[l 0 { k ,jk)lo( k ,jk) + p(jk%( k Jk) + loi k Jk)p{jk)} e ft 3x3 

i=U 

A "*(*) 

Jf(r) = - E "•(MitmiottM + p(«)] € K 3x3 

j=U 
a M fc > 

v7 2 *( r > s ) = - E ™(it)7r(*)7i(fc) G 7l 3x3 

i=i* 

n s (k) 

S k (r) = E ['»0*)pCi*)Ai(fc)] x 4»(*, J*) - ./(;*)*?(*) G ^ 3x3 

j=lk 


sj(r,») = E € ft 3 * 3 

i=U 

A ”*W 

Kt(r) = E 2?>(t,j'*)l™(j*)#0'0^(t)] x - J(Jk)H{k) + H(k)J(j t ) € ft 3 * 3 

J = l* 

= 2[S*(r)]' - E 

i=U 

n a (fc) 

*i(r, >) = E 2^i(fc) [r«O t )p(j*) A; ( fc)] * = 2[S*(r, »)]• € 7J 3 * 3 
i=u 

*{« = E JUk)HWen 3 

J=U 

^*( r > s ) = E pr(^O'fc) - m (i-t)M^,ifc)A;(A;)pOfc)]A^(^) € 7l 3 

j=U 
«»(*) 

= E - m (jk)i J Q ( k )K( k )p(jk)K{ k ) g ft 3 

i=ifc 
nj(fc) _ 

W?(r,s) = E K( k ) m (jk)p(jk)li{ k ) G 7l 3 

i=H 


$ WW*) + JUk)K(i) 


44 



W' 2 (r,s) & E AJWJWJ'W e R 3 

j=l* 

A n ‘ W 

L k (r,s) = ~[l/m(k)} E rn(j k )\ 3 r (k)p(j k )X’{k) € K 3 

j=U 

A n *W . . . . 

T k (r,s) = E [ m Uk)i 3 r { k )p(jk) + J(jk)K( k ))K( k ) € 
j = u 

»«(*) 

j?(r, ») = e [myo#y*)^w - *;<*) jy»)]Ai(i) € r 3 

i=U 

n.(Jfc) , 

= ^[V q { k )]'{Mjk)l J r( k )p(jk) + J(jk)K( k )]K( k ) ZK 1 ( B - 4 ) 

i=U 


Note that 

<?*(r,s) = G fc (s,r) and J k {r,s) = J k (s,r) 

Also define, 


p(fc) t 
F k (r) t 
N k (r) = 

- 7 (*) ± 

S k (r) ± 

K k ( r ) = 


R k (r,s) ± 


Po+Y Pi( s )v( s ) G R 3 

5 = 1 

n m (^) 

f?M + E 

5 = 1 
n m (fc) 

[J,*(r)+ E J‘(v)>lWrsR M 

5=1 

n m (Ar) n m (fe ) n m(^ ) 

J 0 * + E [^"( r ) + {*7i*(0}*M0 + E E «#MM r M s )€ft 3x3 

r= 1 r=l 5=1 

n m (k) 

£f( r ) + E 5 2( r > 5 ) 7 ?( 5 ) G ft 3x3 

5=1 

(^) 

A»+ E *J(r,»)9(«)eR 3<3 

5 = 1 
nm(^) 

fl£(r,s) + E G ( B - 5 ) 

9=1 
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B.2 Modal Mass Matrix 

We have from Eq. (3.25) that the modal mass matrix of the k th body is given by 
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The superscript i = 0,1,2 in denotes the order of dependency of the terms on the 

deformation variables. 







B.3 Expression for a m (k ) 


In this section we derive explicit expressions for the Coriolis and centrifugal acceleration 
term a m (k ). Since 




f 0 l(x,y) \ 

0 ) 


it follows from Eq. 

/ 

£(* + 1 ,*) = 

V 

( 


\ 


(3.14) and Eq. (B.l) that 

o mk + i)Ykhn,k) + [n\k + i)Y<f>(t k+1 ,^ y 

0 <j>(k + l,k) j 

0 [X*(k + 1)]* [A \k + l)]*/(4+i, k) + [Y(k + 1)]* + [A l (k + l)]*J(i fc+1 , k) 
0 0 l(k+\,k) 

0 0 0 


Recalling that the spatial velocity of frame Tk is 


V(k) = 


1 w(Jb) 

V v ( k ) 


where w(k) and v(k) denote the angular and linear velocity respectively of T k we have that 


U j {k) = 

f X j {k) } 


( u{k)\ j {k) ^ 

— 


V V( k ) ) 


\ w(k)y(k) j 


48 


And thus 


/ 


0 


\ 


$'(Jfc + l,*)V m (* + l) = 


tjj(k + l)6„(t k+l ) 


— /(tjt+i, k)u(k + l)^u>(tfc+i) + &(k + l)^u(tfc+i) 

+$w(h+i)i{tk+u k) + d/(k + + 1> k) j 


The vector above has been partitioned so that the term on the top corresponds to modal 
accelerations, the term in the middle to the angular acceleration and the term at the bottom 
to the linear acceleration of the body. Also 


i(k + ljk+i) = &(& + l)l{k + 1, jfc + i) + 6„(j*+i) 


and 


k) = uj{t)c+l)l(tk+l J k) + A„(fc) — + & w (k)l(Oki k) 

— [u)(£ + 1) + ^L(^fc+ i)]^(^Jt+i» k) T A„(&) — &v{dk) + Au{k)l{Oki k) 


where 


Av(fc) 


' A„(fc) ^ 

^ A v (k) J 


H*(k)0(k) 


(B.9) 


Thus 


l{k + 1, k) 


l(k -T 1, tk+\) + ^Jfe+i? k) 

u(k T l)/(& + 1, k) + fi v {tk+i) + ^w(^fc+i)^(^fc+ii k) + A v (k) 
-6 v (d k ) + A„(k)l(O k ,k) 
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Also 


and 



' 0 0 ^ 


H'(k) 

i(o t ,k) 


A(O k ) 

0 


0 

A(O k ) 


\ 

H*(k) 

) 


u(Ok)l(Ok, k) - S v (d k ) 


Thus we have that 


( 


0 


\ 


n*(k) x (k) = 


A(0 k )A w (k)-A(k)6 u (d k ) 


\ u(0 k )A v (k)-u,(k)6 v (d k ) + A w {k)i{O k ,k) -l{O k ,k)u{O k )A„{k) ) 
From Eq. (3.27) and the above expressions it follows that 


/fX d$*(fc+l,Jfc) T/ „ , x dW{k) ... 
a m (&) = — K«(^ + l)d J ~ 4 X(^) — 


dt 


dt 


(B.10) 


a m R{k) 


where 


A(k+l)6Uh+x)+A(k)AUk)-u;(O k )S w (d k ) 


a m R(k) 


u{k + 1) [ut(k + 1 )l(k + 1 ,k) + 2^(t fc+1 )] + [w(4+i) + A{O k )] [v(k) - v(Ot)\ 

+ Jc5(A: + 1) + a)(i* ;+ i)j 3u,(t k +\)l(t k , k — 1) 

~[2A u (k) - 6„(d k )]6 v (d k ) 
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^ £j(k + 1 )u>(k + l)l 0 (k + 1, k) + [u(k + 1) + £(A:)]KA;) - v(Ot)] ) 

- - - _ 


u>{k + l)£up(£k+i) — (^{k)8 u/ {dk) 


+ 


U)(k + 1 ) j^o> ( A: + l)[<5/(£fc+i) — £/(djt)] + 2^ v (ffc + i)| 

+[<L(£fc+i) - L{dk)][v{k) - v(0 ^)] 

+2uj(k + l)^L(£fc+i)A)(£fc+i7 k) 2A w (k)S v (dk) J 


«U(*) 


+ 


2u>(k + 1 )^L(^fc+l)^/(^fc+l) + fiw{tk+l)f>w{tk+l)lo{tki k — 1) + Su>(dk)fiv(dk) J 

~ ” ~ Sw ™~ ~ " ” 


+ 


^ fiu/(tk+l)£>uj(tk+l)fil(tk, k — 1 ) 


(B.ll) 


In the above, denotes the deformation independent part of the Coriolis acceleration, 

while a^ R (it), a^^k) and a^ R (fc) denote the parts whose dependency on the deformation is 
up to first, second and third order respectively. 


B.4 Expression for b m (k ) 

We have from Eq. (3.30) that 

a(h) = dmh), j-(k <n )} v ^ k) = Aj(iW(t) + 4,-(k jk) v(k) 


51 


Since, 


= &(k)l{k y jk) + 6 v (jk) 


it follows that 




a tik) = 


Also from Eq. (3.33) we have that 


Kjk) = 



™(jk)u{jk)u(jk)p(jk) 


Thus, 


b(j k ) + M s (j k )a(j k ) 


u(jk)J(jk)u(jk) + J(jk)^(k)L(jk) + rn(jk)pUk)u(k)[u;(k)l(k,jk) + 26 v (j k ) 


- ptikWfySudk) + L(jk)L(jk)p(jk) + L(jk)Cj(k)p(jk) 

+Cb(k) \ij(k){l(kjk) + p(jk)} + L(jk)p{jk) + 2S v (j k )] | 


(B.12) 
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From Eq. (3.39) we write 


( n-(fc) t a 

U*4 = [b(k) + M,(k)a(k)] £ 

\ m j 


V *t } 

We develop expressions for b*(r), b * and 6* in Eq. (B.13) below. From Eq. (B..12) and 
Eq. (B.13) we have that 

n s (t) 

tjw = e -s(k)\i(k)j(j k )u,(k) - u-(k)\i(k)j(j k )s„(j k ) 

j=U 

-Lu*rk(k)ju 

-sM'KWJUk)^(ik) - [Aj(i)l ■J(h)LUkMk) 

+w-(k){m(j k )p(j k )\i{k)}* [ - + 2«„(j*)] 

+^(iit)[7r(^)]*p(iife)^(ifc)w(fc) 

-m(j fc )a;*(/:) 7 ^(Ar)[ - {l{kj k ) + p(jk)}^(k) + 2^(j fc )] 

+m(j k )hi(k)Y{L{jk)L(jk)p{jk) + L(jk)u(k)p(j k ) + u(k)L{jk)p(jk)} 


n .(*) 



(B.14) 

- w *(^[m(i fc )p(i fc )Ai(fc)] x /(A:,i,V(^) 

(B.15) 

+w*(fc)m(j fc ) 7 ^(A:){/(fc,jfc) + pO’*)M*0 

(B.16) 

+2u;*(fc)(m(jjt)p(i*)A^(fc)] x ^ v (ifc) 

(B.17) 

-2a;*(fc)m(; fc )7r(^v(iifc) 

(B.18) 

-W(*)]V(j*)k&M*) 

(B.19) 
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(B.20) 

+ ^rO’fc)V(ifc)^( fc M fc ) 

(B.21) 

-SM'HikWhMk) 

(B.22) 

+rn{j k )8„(j k )*il{k)p{j k )8 w (j k ) 

(B.23) 

+8UjkTJ(j k )K(k)8Ujk) 

(B.24) 

+Mjk)h}(k)]*6 u ,{jk)u(k)p(j k ) 

(B.25) 

+™(jk)b}(k)]*u(k)6„(j k )p(j k ) 

(B.26) 


Using the modal integrals defined in Section B.l, the above terms can be expressed in the 
following manner: 

-[5.17] + 5.22 = -«•(*) £ T*(s,r)rj(s) 

1 3 = 1 

■t 

-[5.17] + 5.19 + 5.21 = -w*(fc) £ [r 2 fc (s,r) + VUj fc (r,s) + U^’(s,r)]?7(s) 

5.14 + 5.15 = -v m (k)S k (r)w{k) 

5.16 = -w*(jfc)iV fc (rMfc) (B.27) 

5.23 + 5.24 = ^ 5] r 3 fc (4,r,s)77(g)j?(s) 

9 = 1 3=1 

5.20+5.25 = 5.26 

n m (^) 

5.18 + 5.20 + 5.25 + 5.26 = -2 u*(k) £ 5*(s,r)^(s) 

3=1 

Using these, it follows that 

n m (k) 

= -w'(*)[S*(r) + Ar*(r)LO'it)+ E E *?(«.'•. 

9 = 1 3 = 1 

"m(U 

-w*(*) £ [7f(s,r) + T k (s,r) + W l k (r,s) + W k (s,r) + 2F k (s,r)]r)(s) 

j=3 

= -u‘(k) [s*(r) + JV*(r)]“0't) -“■(<••) E <?*(’•.*«(*)+ E E 7 ?(«.'-.«W(«)»W 

j=3 9=1 5=1 
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(B.28) 


where 

<?*( r, s ) = 7f (J, r) + r‘( 3 , r) + .) + W*(j, r) + 2F‘(s, r) (B.29) 

Once again from Eq. (B.12) and Eq. (B.13) we have that 

"»(*) r i 

bt = u(jk)J(jk)u(jk) + J(jk)u{k)8 w (j k ) + m{j k )p(jk)u{k)[u}(k)l(k,j k ) + 2<5 v (j*) 

+w(fc)[w(fe){/(fc,i*) +p(jA:)} + 2<5 v (;'fc) 

= [^O’fc) - m tik)(p(jk)i{k,j k ) + ~l{jk)p{k,jk) + Hjk)i{k>jk))] <*>(&) (B.30) 


-2mO'fc)pO*) + p(i*)]^(iftM fc ) ( B - 31 ) 

( B - 32 ) 

(B.33) 

+Kjk) m Uk)p(jk)8w(jk)u{k) (B.34) 

+u}{k)J(j k )8 w (j k ) (B.35) 

+ ~8 w (j k )J(j k )8„(j k ) (B.36) 

+l{jk)m{jk)L(jk)L(jk)p(jk) (B.37) 

+m(j k )l(j k )L(j k )u(k)p(j k ) (B.38) 

+rn(j k )l(j k )u/(k)8 w (j k )p(j k ) (B.39) 


Once again, using modal integrals, the above terms can be reexpressed in the following 
manner: 


5.30 = w*(k)J(k)ui(k) 


B. 31 = 


n m {k) 

2 E N k (r)r,(r) 


T - 1 


u(k) 
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n m (fc) 

£. 32 + B . 33 + £.34 + £.38 + £.39 = £ A*(r)iJ(r)u;(lb) 


n m (fc) 

£.35 = w(/t) £ Jtf(r)jj(r) (B.40) 


n m (^) n m (fc) 

£.36 + £.37 = ^ ^ £ fc (r,5)^(r)»)(s) 

r=l 5=1 


This results in the following expression 


n m (fc) 

= w(fc)J(/:)u;(fc)+ £ [2N k (r) + /< fc (r)] fi(r)u(k) 


£m(^) w m(&) 

+*(*) E *fM»M + E E R k (r,s)ri(r)ri(s) (B.41) 

r=l r=l 3=1 


Using Eq. (B.12) and Eq. (B.13) it also follows that 


n a (k) 

b k v = 5Z - m Uk)p(jk)v(k)6 w (jk) + m(j k )u(k)[u(k){l(k,j k ) + p(j k )} + 2£ v (j*)] 
i=U 

+m(j k )6 ul (j k )6 u ,(j k )p(j k ) + m(jk)L(jk)^(fc)p(j k ) + m(j k )£j(k)6„(j k )p(j k ) 


"•W 

i=U 

(B.42) 

+m(j k )u(k)u(k){l(kj k ) + p(j k )} 

(B.43) 

+m(j k )L(j k )L(j k )p(j k ) 

(B.44) 

+2m(j k )u(k)6 v (j k ) 

(B.45) 

+m(j k )6 w (j k )u)(k)p(ji t) 

(B.46) 

+m(; J t)w(fc)^(ifc)p(jfc) 

(B.47) 


Using the modal integrals we have that 

£.43 — m(k)uj(k)uj(k)p(k) 

71m (*) n m (k) 

£.44 = m(k) ^ £ L ( r i s )p( r )v( s ) 

r = 1 3=1 
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(B.48) 


W m (&) 


5.42 + 5.45 + 5.46 + 5.47 = 2 w(fc) £ £*(r>K r ) 


r=l 


and thus 


b: = 


n m (fc) n m (fc) n m (fc) 

m(k)<Z(k)i,(k)p(k) + 2d,(k) £ £‘W«r) + m(i) £ £ L( r, «)i}(r)ii(.) 


r=l 


r=l 3=1 


(B.49) 


Putting together Eq. (B.28), Eq. (B.41) and Eq. (B.49) we have that 


+ (!)]«« ) 

-w*(fc)[5f(n m (Ar)) + J,*(n m (fc))]a>(fc) 


M*) = 




m(k)v(k)u>(k)po(k) 


TXT 


+ 


E"” (t| [«*(1,.WW + {SJ(M) + #(!,*)}>?«] 

-*>■(*) ES!*’ [«*(»-(*), «WW + {Sj(n-(*),») + »)}?(')] 


ES,“’ <i( 


?(<•)]• + 


V 


+w(fc)i?f(r)7)(r) 


u)(fc) + 2 E k r) 
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I 


E n m (^) 
7 = 1 


J2 n "( k) T*(q,l,s)ri(q)n(s) 


\ 


g=l 


E"={ k) T 3 k (q,n m (k),s)?}(q)jj(s) 


+ 


E nm{k) (^) 

r=l Z^a-1 


u(k)J k (r, s)u>(k)q(r)r/(s) + jf(r, s)u(k)q(r)q(s) 


+R%(r,s)rj(r)fi(s)\ 


m(k) Er=i k) Es=i k) L(r,s)q{r)q(s) 

" S*i 


/ 


o 


\ 


0 


+ 


(B.50) 


E ^m (*) (fc) m (*) 

7=1 <Lvr = 1 2^5 = 1 




V 


o y 

*?S) 


58 


